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Abstract. Recently, we have proposed that the interaction between relativistic protons resulting from Fermi first order accel- 
eration in the superbubble of a stellar OB association or in other nearby accelerator and ions residing in single stellar winds of 
massive stars could lead to TeV sources without strong counterparts at lower energies. Here we refine this analysis in several 
directions. We study collective wind configurations produced by a number of massive stars, and obtain densities and expansion 
velocities of the stellar wind gas that is to be target of hadronic interactions. We study the expected y-ray emission from these 
regions, considering in an approximate way the effect of cosmic ray modulation. We compute secondary particle production 
(electrons from knock-on interactions and electrons and positrons from charged pion decay), and solve the loss equation with 
ionization, synchrotron, bremsstrahlung, inverse Compton, and expansion losses. We provide examples where configurations 
can produce sources for GLAST satellite, and the MAGIC/HESS/VERITAS telescopes in non-uniform ways, i.e., with or with- 
out the corresponding counterparts. We show that in all cases we studied no EGRET source is expected. Finally, we comment 
on HESS J1303-631 and on Cygnus OB 2 and Westerlund 1 as two associations where this scenario could be tested. 



1. Introduction 

Early-type stellar associations have long been proposed as 
cosmic rays acceleration sites (Cesarsky & Montmerle 1983, 
Manchanda et al. 1996, Bozhokin & Bykov 1994, Parizot et al. 
2004, also Dorman 1999, Romero et al. 1999). For instance, it 
is expected that collective effects of strong stellar winds and su- 
pernova explosions at the core of the associations will produce 
a large-scale shock (the supperbubble region) which will accel- 
erate particles up to energies of hundreds of TeV (e.g., Bykov 
& Fleishman 1992a,b; Bykov 2001). Such relativistic particles, 
if colliding with a dense medium, may produce significant y- 
ray emission, mainly through hadronic interactions. 

O, B, and WR stars, lose a significant fraction of their 
masses in their winds. Indeed, the ultimate result of a stellar 
wind with a high mass-loss rate is to give back gas mass to the 
interstellar medium (ISM). For the sake of this discussion, we 
make here a distinction between the mass that is yet contained 
in single or collective winds of massive stars (in movement) 
and the mass that is free in the ISM (at rest). If star forma- 
tion is on-going, the latter would greatly dominate, since only 
a fraction of the total gas mass contained in the association is 
to end in stars. However, when the star formation is coeval and 
is currently ended, and particularly if one or several supernova 
explosions have pushed away the free gas mass of the region, or 
when the stars under consideration are located in the outskirts 
of a larger association, the mass contained in the innermost 
regions of the winds can exceed that contained in a similarly 



sized area of the ISM. When computing hadronic y-ray lumi- 
nosities, the mass in winds can not be considered negligible in 
these situations. 

Consequently, the winds of a group of massive stars, partic- 
ularly if located close to a cosmic ray acceleration region, may 
act as an appropriate cosmic ray target (e.g., Romero & Torres 
2003). And, because of the wind modulation of the incoming 
cosmic ray flux, what we discuss in more detail below, only 
high energy particles will be able to penetrate into the wind. 
Thus, only high energy photons might be generated copiously 
enough as to be detected. The proposed scenario may then pre- 
dict new potential sources for the new generation of ground- 
based Cerenkov telescopes, which would at the beginning be 
unidentified due of their lack of lower energy counterpart. 1 In 
view of the several unidentified sources already detected by 
HESS in the hundred of GeV - TeV energy regime (Aharonian 
et al. 2005), the aim of this work is to further evaluate this pos- 
sibility. 

2. The gas within a collective wind 

We adopt a similar modelling to that of Canto et al. (2000) 
(see also Chevalier & Clegg 1985, Ozernoy et al. 1997, Stevens 
& Hartwell 2003) to describe the wind of a cluster (or a sub- 
cluster) of stars. This is a hydrodynamical model that does 

1 For an alternative scenario for producing sources in the GeV- TeV 
regime in non-uniform ways see Bosch-Ramon et al. (2005). 



2 



Domingo-Santamaria & Torres: Hadronic processes in winds 



considers the effects of magnetic fields, see next section for 
further discussion on the expected (low) values of the magnetic 
field in the collective wind region. Consider that there are N 
stars in close proximity, uniformly distributed within a radius 
R c , with a number density 



_3N_ 



(D 



Each star has its own mass-loss rate (Mi) and (terminal) wind 
velocity (V,), and average values can be defined as 



V w = 



1 N 

— Y M t . 
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NM W 



(2) 
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All stellar winds are assumed to mix with the surrounding ISM 
and with each other, filling up the intra-cluster volume with a 
hot, shocked, collective stellar wind. A stationary flow in which 
mass and energy provided by stellar winds escape through the 
outer boundary of the cluster is established. For an arbitrary 
distance R from the center of the cluster, mass and energy con- 
servation imply that 



4ttR 2 pV , 
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where p and V are the mean density and velocity of the cluster 
wind flow at position R and h is its specific enthalpy (sum of 
internal energy plus the pressure times the volume), 



h 



7 P 



y-lp 
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with P being the mean pressure of the wind and y being the 
adiabatic index (hereafter y = 5/3 to fix numerical values). 
From the mass conservation equation we obtain, 



nM w 



whereas the ratio of the two conservation equations imply 
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The equation of motion of the flow is 

dV dP 
pV dR=-dR- nM » V > 

which, introducing the adiabatic sound speed c, 

2 P 

c = y— , 
P 

can be written as 
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From the definition of enthalpy and Equation the adiabatic 
sound speed can be expressed as 



J- 
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Using Equation 0, its derivative dp and Equation dl 2l i in dl It 
one obtains 
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which can be integrated and expressed in more convenient di- 
mensionless variables (v = V/V w and r = R/R c ) as follows 



l + ^±Iv 2 



y- 



l 
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with A an integration constant. 

When R > R c , i.e., outside the cluster, by definition n is equal 
to 0, and the mass conservation equation is 



M, 



4-7T . 
T 



, SSO c = —RcnM H ,=47TR^pV, 



(15) 



where the middle equality gives account of the contribution of 
all stars in the association, and M assoc = 2/ M t is the mass-loss 
rate at the outer boundary R c . 

Substituting Equation (I12> and (I15> into the n — realization 
of Equation dl It one obtains 



dR 
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~V 



(y-m 2 -(y+i)v 2 



2{y-\)(V 2 w -V 2 ) 



(16) 



and integrating, the velocity in this outside region is implicitly 
defined from 



v(l-v^) 



2xl/( r -l) 



= Br 



(17) 



with B an integration constant. 

Having constants A and B in Equations 14l > and dl7l i. see be- 
low, the velocity at any distance from the association center can 
be determined by numerically solving its implicit definitions, 
and hence the density is also determined, through Equation Q 

or 02)- 

From Equation Mil , two asymptotic branches can be 
found. When r — > oo, either v — > (asymptotically subsonic 
flow) or v — > 1 (asymptotically supersonic flow) are possible 
solutions. The first one (subsonic) produces the following lim- 
its for the density, the sound speed and the pressure 
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From the latter Equation, if P m is the ISM pressure far from the 
association, the constant B can be obtained as 



(11) B = 



y - l M, 



2y AnP^R 2 



(21) 
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The velocity of the flow at the outer radius r 
from Equation (I17> 

v r=1 (l-v r=1 2 ) 1 /(r-D = B, 
and continuity implies that 

n -(3y+l)/(5y+l) 



1 follows 



(22) 
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Equation d!4t implicitly contains the dependence of v with r 
in the inner region of the collective wind. Its left hand side is 
an ever increasing function. Thus, for the equality to be ful- 
filled for all values of radius (0< r <1), the right hand side of 
the Equation must reach its maximum value at r- 1 . Deriving 
the right hand side of Equation fl!4i . one can find the velocity 
which makes it to be maximum 



( y - 1 
\7+ 1 



1/2 
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Since v grows in the inner region, the maximum velocity is 
reached at r — 1, and from Equation ( 1221 . 
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With the former value of B, and from Equation (12 It . if 
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the subsonic solution is not attainable (continuity of the veloc- 
ity flow is impossible) and the supersonic branch is the only 
physically viable. In this regime, the flow leaves the boundary 
of the cluster R c at the local sound speed v max (equal to 1/2 for 
y = 5/3) and is accelerated until v = 1 for r — > oo. 

Fig.^shows four examples of the supersonic flow (veloc- 
ity and particle density) for a group of stars generating different 
values of M assoc , V w , and R c , as given in Table 1 . The total mass 
contained up to 10 R c is also included in the Table. A typical 
configuration of a group of tens of stars (see Appendix) may 
generate a wind in expansion with a velocity of the order of 
1000 km s _1 and a mass between tenths and a few solar masses 
within a few pc (tens of R c ). We consider below hadronic inter- 
actions with this matter. 

3. Modulation and counterparts 

As it happens in the solar system for cosmic rays with less en- 
ergy than ~ 100 MeV, not all cosmic rays will be able to enter 
into the collective wind of several massive stars. The differ- 
ence between an inactive target, as that provided by matter in 
the ISM, and an active or expanding target, as that provided 
by matter in a single or a collective stellar wind, is given by 
modulation effects. Although wind modulation has only been 
studied in detail for the case of the relatively weak solar wind 



Table 1. Examples of configurations of collective stellar winds. 
The mass is that contained within 10 R c . «o is the central den- 
sity. 



Model 




v w 


R c 


n 


Wind mass 




[M yr-'] 


[km s" 1 ] 


pc 


cnr 3 


[M G ] 


A 


ltr 4 


800 


0.1 


210.0 


0.13 


B 


1(T 4 


800 


0.3 


23.3 


0.39 


C 


5 x 10~ 5 


1000 


0.2 


20.9 


0.11 


D 


2 x 10- 4 


1500 


0.4 


13.9 


0.56 


E 


2 x 1(T 4 


2500 


0.2 


33.5 


0.17 



(e.g. Parker 1958, Jokipii & Parker 1970, Kota & Jokipii 1983, 
Jokipii et al. 1993), and a proper treatment would have to in- 
clude a number of different effects like diffusion, convection, 
particle drifts, energy change, and terminal shock barriers, a 
first approach to determine whether particles can pervade the 
wind is to compute the ratio between the diffusion and convec- 
tion timescales 



e= *d = (3R 2 /D) 
tc (3R/V{R)) 



(28) 



Here D is the diffusion coefficient, R is the position in the wind, 
and V is the wind velocity. Only particles for which e < 1 will 
be able to overcome convection and enter into the wind region 
to produce y-rays through hadronic interactions with matter 
residing there. A similar approach has also been followed by 
White (1985) when computing the synchrotron emission gen- 
erated by relativistic particles accelerated in shocks within the 
wind. In order to obtain an analytic expression for e for a par- 
ticular star we consider that the diffusion coefficient within the 
wind of a particular star is given by (White 1985, Volk and 
Forman 1982, Torres et al. 2004) 

D ~ l -A r c, (29) 

where A r is the mean-free-path for diffusion in the radial direc- 
tion (towards the star). The use of the Bohm parameterization 
seems justified, contrary to what happens in the solar helio- 
sphere, since we expect that in the innermost region of a single 
stellar wind there are many disturbances (relativistic particles, 
acoustic waves, radiatively driven waves, etc.). In the case of 
a collective wind, the collision of individual winds of the par- 
ticular stars forming the association also produce many distur- 
bances. In any event, a change in the diffusion coefficient (say 
a flatter dependence on the energy E) will affect the value of 
the minimum particle energy that protons need to enter into 
the interacting region (see below). We have proven that un- 
less changes in E m [ n are extreme results are not significantly 
affected. 

The mean-free-path for scattering parallel to the magnetic 
field (B) direction is considered to be A\\ ~ 10r g = lOE/eB, 
where r g is the particle gyro-radius and E its energy. In the 
perpendicular direction A is shorter, A ± ~ r g . The mean-free- 
path in the radial direction is then given by A r = A ± 2 sin 2 6 + 
^H 2 cos 2 6» = r g (10 cos 2 0+sin 2 8), where cos~ 2 6 = 1 +(B </) /B r ) 2 . 
Here, the geometry of the magnetic field for a single star is 




Fig. 1. Examples of configurations of collective stellar winds. Main parameters are as in Table 1. 



represented by the magnetic rotator theory (Weber and Davis 
1967; see also White 1985; Lamers and Cassinelli 1999, Ch. 9) 



and 
B 



■--Mir 



(30) 
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where V* is the rotational velocity at the surface of the star, and 
B+ the surface magnetic field. Near the star the magnetic field is 
approximately radial, while it becomes tangential far from the 
star, where A r is dominated by diffusion perpendicular to the 
field lines. This approximation leads — when the distance to the 
star is large compared with that in which the terminal velocity 
is reached, what happens at a few stellar radii — to values of 
magnetic field and diffusion coefficient normally encountered 
in the ISM. 

Using all previous formulae, 
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Equation d32t defines a minimum energy below which the par- 
ticles are convected away from the wind. E^ir) is an increas- 
ing function of r, the limiting value of the previous expression 
being 



E mm (r » /?*) 



3eB i ,V oa R+ I V* 

iMo.i 



UOG, 



R* 

12R C . 



TeV. 



(33) 



Therefore, particles that are not convected in the outer regions 
are able to diffuse up to its base. Note that E mm (r » R+) 
is a linear function of all B* and V*, which is typically 
assumed as V* ~ 0.1 Voo (e.g., Lamers and Cassinelli 1999). 
There is a large uncertainty in these parameters though, which 
is safe to consider about one order of magnitude. The values 
of the magnetic field in the surface of O and B stars is a topic 
of lively discussion. Despite deep searches, only 5 stars were 
found to be magnetic (with sizeable magnetic fields in the range 
of B+ ~ 100 G) (e.g., Henrichs et al. and references therein), 
typical surface magnetic fields of OB stars are then presumably 
smaller. 

In the kind of collective wind we analyzed in Section 2, 
a first estimation of the order of magnitude of the energy scale 
E mm can be obtained as follows. We consider that the collective 
wind behaves as that of a single star having a radius equal to R c , 
and mass-loss rate equal to that of the whole association, i.e., 
A^assoc- The wind velocity at R c , V* is given by Equation (I24i . 
The order of magnitude of the surface magnetic field (i.e., the 
field at R — R c ) is assumed as the value corresponding to the 
normal decay of a single star field located within R c , for which 
a sensitive assumption can be obtained using Equations J30I 
and i3H , <9(10~ 6 ) G. This results, for the whole association, in 

[E min (r » fl*)] assoc ~ 0.8 (=^H TeV. (34) 



\ lyuG AO.lpc/ 

The value of magnetic field is close to that typical of the ISM, 
and is bound to be consider as an average (this kind of mag- 
netic fields magnitude was also used in modelling the uniden- 
tified HEGRA source in Cygnus, see below and Aharonian et 
al. 2005b). In particular, if a given star is close to R c its con- 
tribution to the overall magnetic field near its position will be 
larger, but at the same time, its contribution to the opposite re- 
gion (distant from it 2 R c ) will be negligible. In what follows 
we consider hadronic processes up to 10 - 20 R c , so that a value 
of magnetic field typical of ISM values is expected. We shall 
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consider two realizations of [E mm (r » ^*)] a 
1 TeV. 



100 GeV and 



4. y-rays and secondary electrons from a cosmic 
ray spectrum with a low energy cutoff 

The pion produced y-ray emissivity is obtained from the neu- 
tral pion emissivity as described in detail in the appendix of 
Domingo-Santamaria & Torres 2005. 

4.1. The normalization of the cosmic ray spectrum 

For normalization purposes, we use the expression of the 
energy density that is contained in cosmic rays, u>cr = 
f E N(E) E dE and compare it to the energy contained by cos- 
mic rays in the Earth environment, u>cr,<b(E) = J E N @ (E) E dE, 
where N p(B is the local cosmic ray distribution obtained from 
the measured cosmic ray flux. We assume the Earth-like 
spectrum to be J&(E) = 2.2£' G ^ 5 cirr 2 s~ I sr~ 1 GeV~ I (e.g. 
Aharonian et al. 2001, Dermer 1986), so that u)cr,®(E > 
lGeV) ~ 1.5eVcm -3 . This implicitly defines an enhancement 
factor, g, as a function of energy 



J E N(E)EdE 



(35) 



We assume that N(E) is a power law of the form N(E) = 
K p E~ a . Values of enhancement » 100 at all energies are typi- 
cal of star forming environments (see, e.g., Bykov & Fleishman 
1992a,b; Bykov 2001; Suchkov et al. 1993, Volk et al. 1996, 
Torres et al. 2003, Torres 2004; Domingo-Santamaria & Torres 
2005) and they would ultimately depend on the spectral slope 
of the cosmic ray spectrum and on the power of the accelerator. 
For a fixed slope, harder than that found in the Earth environ- 
ment, the larger the energy, the larger the enhancement, due to 
the steep decline (oc E~ 2J5 ) of the local cosmic ray spectrum. In 
what follows, as an example, we consider enhancements of the 
full cosmic ray spectrum (for energies above 1 GeV) of 1000. 
With such fixed g, the normalization of the cosmic ray spec- 
trum, K p , can be obtained from Equation d35l for every value 
of slope. Note that K p oc g 7 and thus the flux and y-ray lumi- 
nosity, F y and L y , are linearly proportional to the cosmic ray 
enhancement. 



4.2. Emissivities 

We now proceed to the computation of the y-ray emissivity 
produced by a power law spectrum, with a cosmic ray enhance- 
ment of 1000 above 1 GeV. As an example, we have assumed 
a = 2 and 2.3. Fig.|2]shows the results for the whole spectrum 
of cosmic rays (all cosmic rays, i.e., energies above 1 GeV) and 
for cases when the spectrum has a low energy cutoff (100 GeV 
and 1 TeV). For the production of secondary electrons, knock- 
on and pion processes need to be taken into account. Details 
of the computation can be found in Torres 2004. Numerical 
results for the knock on emissivity are shown in Fig. [5] left 
panel for the same spectra we have used before. If the cosmic 
ray spectrum is modulated, the low energy yield of knock-on 




- all cosmic rays 

-cosmic rays above 100 GeV 

■ cosmic rays above 1 TeV 



E - m p [GeV] 

Fig. 2. Contribution of cosmic rays of different energies to the 
hadronic y-ray emissivity. The medium density is normalized 
to 1 cm -3 and the cosmic ray spectrum is proportional to E~ 23 
(black) and E~ zo (grey), with an enhancement of a thousand 
when compared with the Earth-like one above 1 GeV. The nor- 
malization of each spectrum (of each slope) is chosen to respect 
the value of enhancement. In the case of the harder spectrum 
of a — 2.0, we show only the results for the whole cosmic ray 
spectrum, but a similar decrease in emissivity to that of a = 2.3 
can be observed if lower energy cutoffs are imposed. 



electrons gets dramatically reduced. Only when the electron 
energies (E e ) are sufficiently large so that the minimum pro- 
ton energy required to generate them (E™ m ) is larger than the 
modulation threshold, the emissivities obtained with and with- 
out modulation converge. For the charged pion emissivity, as an 
example, we present numerical results for the case of positrons 
in Fig. [3] middle panel. As in the case of the neutral pion pho- 
ton emissivity, a modulated spectrum would produce a much 
smaller electron and positron emissivities at relatively low en- 
ergies. This difference can reach several orders of magnitude 
when comparing with the spectrum obtained when all cosmic 
rays get to interact. 

4.3. Electron energy losses and distribution 



Having the emissivities of secondary electrons we calculate the 
electron distribution solving the diffusion-loss equation. This 
is Tffi ~ IE [b{E)N(E)\ = Q(E), where Q(E) represents all 
the source terms appropriate to the production of electrons 
and positrons with energy E, t(E) stands for the confinement 
timescale, N(E) is the distribution of secondary particles with 
energies in the range E and E + dE per unit volume, and 
b(E) = - (dE/dt) is the rate of energy loss. The energy losses 
considered are those produced by ionization, inverse Compton 
scattering, bremsstrahlung, synchrotron radiation, and the ex- 
pansion of the medium (for a compilation of the relevant for- 
mulae see the appendix in Torres 2004). A number of uncer- 
tain parameters enters into the computation of these losses. 
Most notably, these parameters are the medium density n af- 
fecting bremsstrahlung and ionization losses, the magnetic field 
B affecting synchrotron losses, the photon target field affect- 
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all cosmic rays 

cosmic rays above 1 00 GeV 

cosmic rays above 1 TeV 



- all cosmic rays 
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■ cosmic rays above 1 TeV 
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Fig. 3. The effect of a modulated cosmic ray spectrum (the same as in Fig. [2] with a = 2.3 and 2.0) over the electron knock-on 
emissivity (left) and the positron emissivity (middle) of a medium with density n — 1 cm -3 . Right: Different losses for 
assumed parameters: Curves 1 correspond to ionization losses for n = 100 and 20 cm 3 . Curves 2 correspond to synchrotron 
losses for B — 50 and 200 fiG (see text). Curves 3 correspond to inverse Compton losses for a photon energy density of 20 and 
100 eV cirT 3 . Curves 4 correspond to Bremsstrahlung losses for n — 100 and 20 cirT 3 . Finally, curve 5 correspond to adiabatic 
losses having a ratio V(km s _1 )//?(pc)=300 (e.g., a wind velocity of 1500 km s and a size relevant for the escape of the electrons 
of 5 pc). 



ing inverse Compton losses, and the velocity of the expanding 
medium and the size relevant for escape, affecting the expan- 
sion losses. 

Fig. |3 right panel shows the rate of energy losses for a 
range of parameters. We show results both n = 100 and 20 
cm" 3 and we allow magnetic fields to reach up to 200 fiG in 
the collective wind region. The latter is made in order to fic- 
titiously enhance -on purpose- the synchrotron losses and to 
better shown the dominance of the other losses mechanisms 
over them. Inverse Compton losses are computed using two 
different normalization for the energy density when the pho- 
ton target is considered to be a blackbody distribution with 
T e ff =50000 K. The expansion losses have the form 




where V is the collective wind of the region, and R its relevant 
size. In Fig. [3] right panel, to be conservative, we have chosen a 
ratio V(km s~')/,R(pc)=300 (e.g., a -single or collective- wind 
velocity of 1500 km s _I and a size relevant for the escape of the 
electrons of 5 pc). Fig. [3] right panel shows that the expansion 
dominates the electron losses b(E), as well as the confinement 
timescale t(E), throughout a wide range of energies. In Fig.0] 
we show an example of the resulting secondary electron distri- 
bution obtained by numerically solving loss equation. Only at 
high energies the effect of the cutoff is unnoticeable, whereas 
it greatly affects the production of secondaries (up to several 
orders of magnitude) below 10 GeV. 

As a benchmark, we note that the radio emission from the 
electron population of Fig. 0] is well below the upper limits 
imposed with VLA at the location of the Cygnus unidenti- 
fied HEGRA source (see below for a more detailed discussion, 



10' I i 




(E^-mJ [GeV] 



Fig. 4. Secondary electron distribution obtained by numerically 
solving the loss equation. The primary cosmic ray spectrum is 
the same as in Fig. |3 with a = 2.3, and results are shown for 
different low energy cutoffs. The average density is assumed as 
n — 20 cirT 3 , the magnetic field is assumed as 50 ^G, and the 
size relevant for escaping the modulating region is 5 pc. The 
inset shows the total energy loss rate b(E). 

< 200 mJy at 1.49 GHz, Butt et al. 2003). Even assuming a 
rather high magnetic field of 50 fiG as in Fig.0] what is in the 
case of the HEGRA source discarded by X-ray and radio ob- 
servations, we obtain ~ 50 mJy at the quoted frequency for 
the non-modulated cosmic ray population. A smaller magnetic 
field, as is probably to be found in the outer regions, or a mod- 
ulated production of secondaries will diminish this estimation. 
We verify too that the limit is respected even when considering 
that the primary electron population (particularly at low ener- 
gies) is at the same level than the secondary electron distribu- 
tion. 
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4.4. Total y-ray flux from a modulated environment 

To compute y-ray fluxes in a concrete example, and following 
Section 2, we consider ~2 M Q of target mass being modulated 
within ~1 pc. The average density is ~ 25 cm 4 . This amount 
of mass is typical of the configurations studied in Section 2 
within the innermost 20 R c ~ 2 - 8 pc. To fix numerical values, 
we consider that the group of stars is at a Galactic distance 
of 2 kpc. Using the computations of secondary electrons and 
their distribution, the left panel of Fig.|3]shows the differential 
(hadronic and leptonic) y-ray flux when the proton spectrum 
has a slope of 2.3 and 2.0. In the latter case, to simplify, we 
show only the pion decay contribution which dominates at high 
energies, as is produced by the whole cosmic ray spectrum. 

The differential photon flux is given by F y (E y ) = 
[V/4nD 2 ]Q y (E y ) = [M/m p 4nD 2 ][Q y (E y )/n], where V and D 
are the volume and distance to the source, and M the target 
mass. In those examples where the volume, distance, and/or 
the medium density are such that the differential flux and the 
integral flux obtained from it above 100 MeV with the full 
cosmic ray spectrum is larger than instrumental sensitivity, a 
modulated spectrum with a 100 GeV or a 1 TeV energy thresh- 
old might not produce a detectable source in this energy range. 
However, the flux will be essentially unaffected at higher en- 
ergy (see the high energy end of Fig. s |21 El and@}. The left 
panel of Fig. [5] shows that wind modulation can imply that a 
source may be detectable for the ground-based Cerenkov tele- 
scopes without being even close to be detected by instruments 
in the 100 MeV - 10 GeV regime (like the late EGRET or the 
forthcoming GLAST). To have further insight on this, the right 
panel of Fig. [5] presents the integral flux of y-rays as a func- 
tion of energy, together with the sensitivity of ground-based 
and space-based y-ray telescopes. The sensitivity curves shown 
are for point-like sources, it is expected that extended emission 
would require about a factor of 2 more flux to reach the same 
level of detectability. Table 3 summarizes these results. From 
Table 3 and Fig.|5]it is possible to see that there are plenty of 
different scenarios (possibly relevant parameters are distance, 
enhancement, degree of modulation of the cosmic ray spec- 
trum, and slope) for which sources that shine enough for detec- 
tion at the GLAST domain, may not do so in the IACTs energy 
range, and viceversa. 

5. Opacity to y-ray escape 

The opacity to pair production of the y-rays in the UV stel- 
lar photon field of an association can be computed as (Reimer 
2003, Torres et al. 2004) 

r(£ r ) = 2 I I N i (E ir )(r e -AE*,E y )dE ir dn, (37) 
,•=1 J JRci 

where is the energy of the soft photons, E y is the energy 
of the y-ray, R c j is the place where the photon was created 
with respect to the position of the star number i, r, is the dis- 
tance measured from star number i, and cr e - e +(E+, E 7 ) is the 
cross section for yy pair production (Cox 1999, p.214). Note 
that the lower limit of the integral on e in the expression for 
the opacity is determined from the condition that the center of 



Table 2. Examples of results for detection in different tele- 
scopes when the configuration of collective stellar winds gener- 
ates a target of about 2 M , it is located at 2 kpc, and it is bom- 
barded with a cosmic ray spectrum having an spectral slope 
a = 2.3 and 2.0 enhanced a factor of 10 3 above 1 GeV. The 
full cosmic ray spectrum and different modulated cases, at 100 
GeV and 1 TeV, are shown. Except in the case of EGRET, when 
a = 2.3 sensitivities are barely above the expected fluxes (see 
Fig. 



Telescope 


All 


100 GeV 


1 TeV 


a = 2.3 


EGRET (E > 100 MeV) 


X 


X 


X 


GLAST (E > 100 MeV) 


V 


X 


X 


MAGIC (E > 50 GeV) 


X 


X 


X 


HESS/VERITAS (E > 100 GeV) 


X 


X 


X 


a = 2.0 


EGRET (E > 100 MeV) 


X 


X 


X 


GLAST (E > 100 MeV) 


V 


X 


X 


MAGIC (E > 300 GeV) 


V 


V 


V 


HESS/VERITAS (E > 200 GeV) 


V 


V 


V 



mass energy of the two colliding photons should be such that 
(1 - (mc 2 ) 2 /(e E y )) > 0. The stellar photon distribution of star 
number i at a position r, from the star is that of a blackbody 
peaking at the star effective temperature (T e fff) and diluted by 
the distance factor, V,(£*) = (ttB(£ , *))/(/i£ , *c) • (Rl.)/(r 2 ), 
where h is the Planck constant, R+. is the star number i radius, 
and B(£+) = (2E+ 3 )/((hc) 2 ) ■ - l) -1 . 

The size of the emission region (where the target gas mass 
is located) that we have considered as an example in the pre- 
vious sections is of the order of 1 pc. The typical radius of a 
massive star is 10-20 solar radius ~ 5 xl0~ 7 pc, so that the 
most likely creation sites for photons will be far from individ- 
ual stars. In Fig.|5](right panel) we show the value of r(E y ) for 
different photon creation sites distant from a 04V-star 10, 100, 
and 1000 R@, with fi* = 12R @ and 7/ eff = 47400 K. Unless a 
photon is created hovering onto the star, well within 1000 
y-ray opacities are very low and can be safely neglected. This 
is still true for associations in which the number of stars is of 
some tens. Consider for instance a group of 30 such stars within 
a region of 0.5 pc (the central core of an association). The stel- 
lar density is given by Equation and the number of stars 
within a circle of radius R progresses as N — N(R/R C ) 3 - Fig. 
H] (right panel) shows the collective contribution to the opac- 
ity obtained from Equation (I37> in this configuration, it is also 
very low since the large majority of the photons are produced 
far from individual stars. This, however, is not the case if one 
considers the collective effect of a much larger association like 
Cygnus OB 2, particularly at its central region (Reimer 2003). 
Reimer demonstrated that even when a subgroup of stars like 
the ones considered here is separated from a super cluster like 
Cygnus OB2 by about 10 pc, the influence of the latter pro- 
duces an opacity about one order of magnitude larger than that 
produced by the local stars. Even in this case, Fig. [5] (right 
panel) shows that this opacity is not enough to preclude es- 
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Fig. 5. Differential (left) and integral (right) fluxes of y-rays emitted in a non-modulated and a modulated environment. The bump 
at very low energies in the left panel is produced because we show leptonic emission coming only from secondary electrons. 
Above ~ 70 MeV the emission is dominated by neutral pion decay. Also shown are the EGRET, GLAST, MAGIC, and HESS 
sensitivities. Note that a source can be detectable by IACTs, and not by GLAST, or viceversa, depending on the slope of the 
cosmic ray spectrum and degree of modulation. Right: Opacities to yy pair production in the soft photon field of an 04V-star at 
10, 100, and 1000 R*, and in the collective photon field of an association with 30 stars distributed uniformly over a sphere of 0.5 
pc. The closest star to the creation point is assumed to be at 0.16 pc, and the rest are placed following the average stellar density 
as follows: 1 additional star within 0.1, 2 within 0.25, 4 within 0.32, 8 within 0.40 and 14 within 0.5 pc. 



cape from the region of the local enhancement of stellar den- 
sity. Such is not the case at the very center of Cygnus OB 2. 

6. Candidates: Cygnus OB 2 and Westerlund 1 

Extensive studies of Galactic and LMC/SMC star clusters and 
OB associations suggest that star formation occurs almost in- 
stantaneously (Massey et al. 1995, Leitherer 1999). Typical age 
spreads are about 2 Myr or less. This is short in comparison 
with stellar evolutionary timescales, except for the most mas- 
sive stars, and so a coeval star formation seems appropriate. In 
addition, we are particularly interested in the case in which at 
the position of the stars that are being illuminated by cosmic 
rays there is no current star formation. If it were, the amount 
of gas and molecular material in the ISM would be larger than 
that contained in the winds, and would make the latter a sub- 
dominant contribution in the generation of the total TeV flux. 
In particular, should collision between cosmic rays and nuclei 
in the ISM be dominant, there would be no modulation. This 
fact and the opacity to y-ray escape that is found in the center 
of a very massive cluster points to the best scenario: a sub- 
group of stars located in the outskirts of an association, close 
to an accelerator region, perhaps a SNR, or the association su- 
perbubble itself. The case of Cygnus OB 2 and the TeV source 
detected by HEGRA, TeV J2032+4131 (Aharonian et al. 2002, 
2005b), seems to be a possible realization of this scenario. This 
TeV source, for which no counterparts at lower energies are 
presently identified (Butt et al. 2003), constant during the three 
years of data collection, extended (5.6+1.7 arcmin, ~ 2.7 pc at 
1.7 kpc), separate in about 10 pc (at 1.7 kpc) from the core of 
the association, and coincident with a significant enhancement 
of the star number density (see Fig. 1 of Butt et al. 2003) might 
be suggestive of the scenario outlined in the previous sections. 



TeV J2032+4131 presents an integral flux of about 3% of that 
of the Crab [F y {E 7 > ITeV) = 4.5(+l. 3) x 10~ 13 photons cirT 2 
s _1 ], and a y-ray spectrum F 7 {E y ) = B(£' y /TeV) _r photons 
cm" 2 s" 1 TeV" 1 , where B = 4.7 (±2.1 stat + 1.3 sys ) x 10~ 13 and 
F = 1.9(±0.3 stat ± 0.3 sys ). Within the kind of models studied in 
this work, it would be possible to explain, apart from consis- 
tency with the flux level, spectrum, and variability, why there 
is no detectable TeV source at the central core of Cygnus OB 2 
(large opacity to photon escape), and why there is no EGRET 
(and will not be a GLAST source) or other significant radio or 
X-ray diffuse emission at the position of the HEGRA detection 
(modulation of cosmic rays). 

However, it is worth noticing that the only stars quoted by 
Butt et al. at the position of TeV J2032+4131 are 10 O and 10 
B stars, and not even all of them are within the contours of the 
source. The quoted stars are similar, late O and early B, and 
together produce a mass-loss rate of about 1 x 10~ 5 M G yr . 
This mass-loss rate is low and produces not too dense a collec- 
tive wind, according to the simulations of Section 2. In addi- 
tion, typical distances between these stars are of the order of a 
parsec, so that a central core radius is not well defined. Thus, 
if we are not missing significant stars in this region, something 
we might very well be judging on newest and deeper Chandra 
observations of the region (Butt et al, and Reimer et al. both in 
preparation) as well as on the expected star density, a better ap- 
proach to study the possible contribution of stellar winds to this 
source is just to add that of individual stars (Torres et al. 2004). 
If only the currenlty known stars exist, the produced flux is too 
low to produce the source unless a larger enhancement or a low 
ISM density are invoked. A simple hadronic scenario where 
an enhanced spectrum interacts with more abundant ISM nu- 
clei can not be discarded at this point, since there is no need 
to use modulation to explain why there is no source detected 
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by EGRET (the predicted flux in this model is below EGRET 
sensitivity). However, this model can be tested with GLAST 
observations (see the example of Fig.|5Jl. Also, MAGIC obser- 
vations of this region are crucial. We have yet to wait for further 
data to reach a definitive answer. 

Knodlseder (2003) has compiled a sample of massive 
young star clusters (Cygnus clones) that have been observed 
in the Galaxy (see his table 2). Westerlund 1 (Wd 1) is a 
prominent member. The usually adopted distance to the lat- 
ter is D ~ 1.1 ± 0.4 kpc (Piatti et al. 1998), although Clark 
et al. (2005) favor a distance between 2-A kpc. The known 
zoo of massive stars, clearly a lower limit in each category, 
includes 7 WN, 6 WC, 5 Early transition stars (like LBV, or 
sgB[e]) and more than 25 OB (Clark and Negueruela 2002, 
Clark et al. 2005). A luminous blue variable (W243) appar- 
ently undergoing an eruption event was also found (Clark & 
Negueruela 2004). W243's mass-loss rate alone can be as high 
as 3 - 6 x 10~ 4 M yr~' . The stellar population of Wd 1 appears 
to be consistent with an age of the order 4-8 Myr if the cluster 
is coeval, and with a lower limit for a total mass of about a few 
thousand M , most likely to be around a few xlO 5 M (Clark 
& Negueruela 2002, Clark et al. 2005). Then, Wd 1 is likely to 
be one of the most massive young clusters in the Local Group, 
and in contrast to Cygnus OB2, it is much more compact (a 
radius of about 0.6 pc) -and thus a better target for pointing 
instruments. A few stars a bit far from the center (so that high 
opacities are avoided) subject to illumination by a cosmic ray 
population with a hard slope would make Wd 1 a y-ray source. 
The HESS observatory has covered the position of Wd 1 when 
doing the galactic plane scan, although with just a few hours 
of observation time (Aharonian et al. 2005). We suggest that 
pointed observations to Wd 1 might result in its detection and 
that a model like the one presented here can be tested jointly 
by GLAST and HESS. 

Finally, special interest on this model may appear in re- 
gards to the serendipitous discovered HESS J 1303-631 TeV y- 
ray source (Aharonian et al. 2005c), extended and still uniden- 
tified, which is in spatial coincidence with the OB association 
Cen OB6 (see their Fig. 7), that contain al least 20 known O 
stars and even 1 WR star. 

7. Concluding Remarks 

We have studied collective wind configurations produced by a 
number of massive stars, and obtained densities and expansion 
velocities of the stellar wind gas that is to be target for hadronic 
interactions in several examples. We have computed secondary 
particle production, electrons and positrons from charged pion 
decay, electrons from knock-on interactions, and solve the ap- 
propriate diffusion-loss equation with ionization, synchrotron, 
bremsstrahlung, inverse Compton, and expansion losses to ob- 
tain expected y-ray emission from these regions, including 
in an approximate way the effect of cosmic ray modulation. 
Examples where different stellar configurations can produce 
sources for GLAST satellite, and the MAGIC/HESS/VERITAS 
telescopes in non-uniform ways, i.e., with or without the cor- 
responding counterparts were shown. Finally, we have com- 
mented on Cygnus OB 2 and Westerlund 1 as two associations 



Table 3. Wind model parameters of WR, O, and B stars. 



Stellar 


log[M ] 




type 


M yr" 1 


[km s" 1 ] 


WNL 


-4.2 


1650 


WNE 


-4.5 


1900 


WC6-9 


-4.4 


1800 


WC4-5 


-4.7 


2800 


WO 


-5.0 


3500 


03 


-5.2 


3190 


04 


-5.4 


2950 


04.5 


-5.5 


2900 


05 


-5.6 


2875 


05.5 


-5.7 


1960 


06 


-5.8 


2570 


06.5 


-5.9 


2455 


07 


-6.0 


2295 


07.5 


-6.2 


1975 


08 


-6.3 


1755 


08.5 


-6.5 


1970 


09 


-6.7 


1500 


09.5 


-6.8 


1500 


BO 


-7.0 


1000 


B0.5 


-7.2 


500 


Bl 


-7.7 


500 


Tl 1 ^ 

r> 1 . j 


-o.Z 


ouu 


B2 


-8.6 


500 


B3 


-9.5 


500 


B5 


-10.0 


500 


B7 


-10.9 


500 


B8 


-11.4 


500 


B9 


-12.0 


500 



where this scenario could be tested. In the latter case, we have 
proposed it to be targeted by HESS. 

Appendix: Mass loss rates and terminal velocities 
of individual stars 

To estimate how many stars, and of what kind, are needed to 
get certain averages for the association parameters, we adopt 
the theoretical wind model C of Leitherer, Robert & Drissen 
(1992), where it is discussed in more detail. This model, re- 
ferred to as Model THEOR in Leitherer et al. (1999), is the 
standard model used in the synthesis program STARBURST99 
to compute spectrophotometric and statistics of starburst re- 
gions. 2 We remark that mass-loss rates and terminal veloci- 
ties are indeed uncertain (typically about ~ 30%). Varying the 
mass-loss rate and terminal velocity to those corresponding to 
other models will have an impact in any y-ray luminosity com- 
putation. Alternative parameterizations to the ones below for 
M* and V m can be found, for instance, in the work by Vink. et 
al. (2000). 

For WR stars 3 , mass-loss rates are based on observations 
by van der Hucht et al. (1986) and Prinja et al. (1990). Average 

2 http://www.stsci.edu/science/starburst99 

3 The first evolutionary phase of a WR star is the nitrogen-line stage 
(WNL), which begins when CNO processed material is exposed at the 
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terminal velocities have been taken from Prinja (1990). For O 
and B stars, we shall use the parameterizations of M (M© yr _1 ) 
and Voo (km s _1 ) in terms of stellar parameters (Leitehrer et 
al. 1992, 1999), that is also part of Model THEOR. A multi- 
dimensional fit to M (M Q yr _1 ) and Voo (km s~') as a function 
of stellar luminosity L, mass M, effective temperature T e g, and 
metallicity Z, gives 

logMtMoyr 1 ) = -24.06 + 2.45 log[L(L )] 

-l.lOlog[M(M )] + 1.31 logfTdK*)] 
+O.8Olog[Z(Z )], (38) 

y (X ,(kms- 1 ) = 1.23 - 0.30 log[L(L Q )]+ 0.55 log[M(M Q )] 

+0.641og[r eff (^)] +0.13 log[Z(Z )]. (39) 

In order to obtain the mass-loss rates and terminal velocities, 
we then need to apply a calibration between the spectral sub- 
types of the O and B stars and their stellar parameters. For early 
B and O stars we adopt the work by Vacca et al. (1996), see 
their Tables 5-7. Columns 2 and 3 of Table 1 show the stellar 
wind parameters resulting from the previous calibration. Since 
in general, many of the stars in a young association will pertain 
to luminosity class V, we give the results for this class only, 
although combining this Section with the work by Vacca et al. 
(1996), values for other classes could be obtained as needed. 
We adopt the values of masses as derived using evolutionary 
codes, since as discussed by Vacca et al. (1996), the cause of 
the mass discrepancy between results from the former codes 
and spectroscopic analysis seems to be that the models used 
in the spectroscopic studies do not properly take into account 
the effects of the wind extension, mass outflow, and velocity 
fields, and are thus less reliable. For B stars of late type, we 
complement the calibration by Vacca et al. (1996) with that of 
Humphreys & Mc Elroy (1984), see their Table 2, and Schmidt- 
Kaler (1982), which give the bolometric corrections, {BC), and 
effective temperatures. The total luminosity is then computed 
using the equation 

log(L/Lo) = -0.4(M, + BC - M bol0 ), (40) 

where Mb o i = 4.75 mag (Vacca et al. 1996). The masses and 
visual magnitudes M v of B stars are obtained following the 
method described in Appendix A of Knodlseder (2000). The 
star radius of each of the stars can be obtained by the usual re- 
lationship R\ = LI(4naT A eS ), where cr = 5.67 x 10~ 5 erg s" 1 
cm" 2 KT 4 . 

It is to be noted that this mass-loss parameterization may 
yield to an under-estimation for certain stars, and an over- 
estimation of the terminal velocities, like for example in the 
case of HD93129A (Benaglia & Koribalski 2004). We re- 
mark that almost 400 Galactic O stars have recently been com- 
piled in the new 'Galactic O Star Catalogue' (GOS) by Mafz- 
Apellaniz & Walborn (2002), but there are only about a dozen 

stellar surface. The second stage is the WNE, during which no surface 
hydrogen is detectable. After the helium envelope is shed, the WC 
stage occurs, in which strong carbon lines can be seen, and finally, WO 
stars are formed, with high surface oxygen abundances. Further sub- 
classifications are made according to line strength ratios (e.g., Smith 
& Maeder 1991). 



stars known with spectral types 03.5 or earlier, thus limiting 
the statistical knowledge for comparison. In any case, this pa- 
rameterization is enough to show that a group of several tens of 
early stars can easily reach a M assoc ~ (9(1 -5 - 10~ 4 ) M yr . 
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